// Position.h: interface for the Position and PMetric classes
//
//////////////////////////////////////////////////////////////////////

#if !defined(AFX_POSITION_H__B3C5CB6A_4910_11D1_939A_006097737B41__INCLUDED_)
#define AFX_POSITION_H__B3C5CB6A_4910_11D1_939A_006097737B41__INCLUDED_


#include <cmath>
#include <iomanip>
#include <iostream>
using namespace std;
#include <stdlib.h>


class Position;

class PMetric
{
	friend Position operator*(const PMetric  &M, const Position &Y);
	friend Position operator*(const PMetric  &M, const double   *y);	
	friend PMetric  PMatrix (const Position &X0, const Position &X1, const Position &X2);
	friend PMetric  TPMatrix(const Position &X0, const Position &X1, const Position &X2);
	
public:
	PMetric()	{}
	PMetric(const PMetric &M);
	void operator=(const PMetric &M);



private:
	double e[3][3];

	void Copy(const PMetric &M);
};

//////////////////////////////////////////////////////////////////////
// inline functions
//////////////////////////////////////////////////////////////////////

inline void PMetric::Copy(const PMetric &M)
{
int n,m;
for(n=0;n<3;n++)
	for(m=0;m<3;m++)
		e[n][m] = M.e[n][m];
}

//////////////////////////////////////////////////////////////////////
// inline constructors
//////////////////////////////////////////////////////////////////////

inline      PMetric::PMetric  (const PMetric &M)		
{	Copy(M);	}

//////////////////////////////////////////////////////////////////////
// inline overloaded operators (class members)
//////////////////////////////////////////////////////////////////////

inline void PMetric::operator=(const PMetric &M)	
{	Copy(M);	}



//////////////////////////////////////////////////////////////////////
//
// interface for the Position class.
//
//////////////////////////////////////////////////////////////////////


class Position  
{
	friend class PMetric;  // so PMetric is a friend of Position
                         // and therefore PMetric can access all private
                         // data of Position

	friend Position operator+ (const Position &X, const Position &Y);
	friend Position operator- (const Position &X, const Position &Y);
	friend double   operator* (const Position &X, const Position &Y);
	friend Position operator* (const Position &X, const double   &y);
	friend Position operator* (const double   &x, const Position &Y);
	friend Position operator* (const PMetric  &M, const Position &Y);
	friend Position operator* (const PMetric  &M, const double   *y);	
	friend Position operator/ (const Position &X, const double   &y);
	friend Position operator^ (const Position &X, const Position &Y);
	friend bool     operator==(const Position &X, const Position &Y);
	friend bool		operator!=(const Position &X, const Position &Y);
	friend ostream& operator<<(ostream &os, const Position &X);
	friend PMetric  PMatrix (const Position &X0, const Position &X1, const Position &X2);
	friend PMetric  TPMatrix(const Position &X0, const Position &X1, const Position &X2);

public:
	Position(double x, double y, double z);
	Position(){}
	Position(const Position &X);
   ~Position(){}
	void Normalise(void);
	void Normalise(double Sign);
	void Zero(void);

   // used in class Atom when non-Z-matrix description is used

   void setXYZ(const double x, const double y, const double z) {xyz[0]=x; xyz[1]=y; xyz[2]=z;}
   void setX(const double x) {xyz[0]=x;}
   void setY(const double y) {xyz[1]=y;}
   void setZ(const double z) {xyz[2]=z;}
   double x() const {return xyz[0];};   // used in molecule::TransformDerivatives
   double y() const {return xyz[1];};
   double z() const {return xyz[2];};


	Position operator=( const Position &X);
	void     operator+=(const Position &X);
	void     operator-=(const Position &X);
	void     operator*=(const double   &x);
	void     operator/=(const double   &x);

private:
	double xyz[3];	

};

//////////////////////////////////////////////////////////////////////
// inline constructors
//////////////////////////////////////////////////////////////////////

inline Position::Position(double x, double y, double z)
{
xyz[0] = x;	xyz[1] = y;	xyz[2] = z;
}

inline Position::Position(const Position &X)
{
for(int i=0;i<3;i++)	xyz[i] = X.xyz[i];
}

//////////////////////////////////////////////////////////////////////
// inline functions
//////////////////////////////////////////////////////////////////////

inline void Position::Normalise(void)
{
int    i;
double Norm=0.0;

for(i=0;i<3;i++)	Norm += xyz[i]*xyz[i];
Norm = 1.0/sqrt(Norm);
for(i=0;i<3;i++)	xyz[i] *= Norm;
}

inline void Position::Normalise(double Sign)
{
Normalise();
(*this) *= Sign;
}

inline void Position::Zero(void)
{
for(int i=0;i<3;i++)	xyz[i] = 0.0;
}

//////////////////////////////////////////////////////////////////////
// inline overloaded operators (class members)
//////////////////////////////////////////////////////////////////////

inline Position Position::operator=(const Position &X)
{
for(int i=0;i<3;i++)	xyz[i] = X.xyz[i];
return (*this);
}

inline void Position::operator+=(const Position &X)
{
for(int i=0;i<3;i++)	xyz[i] += X.xyz[i];
}

inline void Position::operator-=(const Position &X)
{
for(int i=0;i<3;i++)	xyz[i] -= X.xyz[i];
}

inline void Position::operator*=(const double   &x)
{
for(int i=0;i<3;i++)	xyz[i] *=x;;
}

inline void Position::operator/=(const double   &x)
{
for(int i=0;i<3;i++)	xyz[i] /=x;;
}

//////////////////////////////////////////////////////////////////////
// inline overloaded operators (friends)
//////////////////////////////////////////////////////////////////////

inline Position operator+(const Position &X, const Position &Y)
{
return (Position(X.xyz[0]+Y.xyz[0],X.xyz[1]+Y.xyz[1],X.xyz[2]+Y.xyz[2]));
}

inline Position operator-(const Position &X, const Position &Y)
{
return (Position(X.xyz[0]-Y.xyz[0],X.xyz[1]-Y.xyz[1],X.xyz[2]-Y.xyz[2]));
}

inline double   operator*(const Position &X, const Position &Y)
{
return (X.xyz[0]*Y.xyz[0] + X.xyz[1]*Y.xyz[1] + X.xyz[2]*Y.xyz[2]);
}

inline Position operator*(const Position &X, const double   &y)
{
return (Position(X.xyz[0]*y,X.xyz[1]*y,X.xyz[2]*y));
}

inline Position operator*(const double   &x, const Position &Y)
{
return (Position(x*Y.xyz[0],x*Y.xyz[1],x*Y.xyz[2]));
}

inline Position operator*(const PMetric &M, const Position &Y)
{
double y0=0.0,y1=0.0,y2=0.0;

for(int n=0;n<3;n++)
	{
	y0 += M.e[0][n]*Y.xyz[n];
	y1 += M.e[1][n]*Y.xyz[n];
	y2 += M.e[2][n]*Y.xyz[n];
	}

return Position(y0,y1,y2);
}

inline Position operator*(const PMetric &M, const double *y)
{
double y0=0.0,y1=0.0,y2=0.0;

for(int n=0;n<3;n++)
	{
	y0 += M.e[0][n]*y[n+1];
	y1 += M.e[1][n]*y[n+1];
	y2 += M.e[2][n]*y[n+1];
	}

return Position(y0,y1,y2);
}

inline Position operator/(const Position &X, const double   &y)
{
return (Position(X.xyz[0]/y,X.xyz[1]/y,X.xyz[2]/y));
}

inline Position operator^(const Position &X, const Position &Y)
{
return (Position(X.xyz[1]*Y.xyz[2] - X.xyz[2]*Y.xyz[1],
				 X.xyz[2]*Y.xyz[0] - X.xyz[0]*Y.xyz[2],
				 X.xyz[0]*Y.xyz[1] - X.xyz[1]*Y.xyz[0]));
}

 inline ostream& operator<<(ostream &os, const Position &X)
{
  os.precision(5);
  os << setw(10) << X.xyz[0] << setw(10) << X.xyz[1] << setw(10) << X.xyz[2];

  return (os);
}

inline bool operator==(const Position &X, const Position &Y)
{
for(int i=0;i<3;i++)
	if(fabs(X.xyz[i]-Y.xyz[i])/(fabs(X.xyz[i])+fabs(Y.xyz[i])+1.0)>1e-8) return (false);

return (true);
}

inline bool operator!=(const Position &X, const Position &Y)
{
return (!(X==Y));
}

//////////////////////////////////////////////////////////////////////
// inline functions (friends)
//////////////////////////////////////////////////////////////////////

inline PMetric PMatrix(const Position &X0, const Position &X1, const Position &X2)
{
PMetric M;

for(int n=0;n<3;n++)
	{
	M.e[0][n] = X0.xyz[n];
	M.e[1][n] = X1.xyz[n];
	M.e[2][n] = X2.xyz[n];
	}

return M;	
}

inline PMetric TPMatrix(const Position &X0, const Position &X1, const Position &X2)
{
PMetric M;

for(int n=0;n<3;n++)
	{
	M.e[n][0] = X0.xyz[n];
	M.e[n][1] = X1.xyz[n];
	M.e[n][2] = X2.xyz[n];
	}

return M;
}
#endif // !defined(AFX_POSITION_H__B3C5CB6A_4910_11D1_939A_006097737B41__INCLUDED_)
